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A study is made of methods of solution of a pair of 
nonlinear, coupled, third-order differential equations in 
two dependent and two independent variables. The equations 
have an incomplete initial value set. Nonlinearity is a 
result of the interdependence of solutions and coefficients. 
These equations are found to be of interest in the study of 


compressible boundary layer theory. 


The development of methods based on a technique of 
representing the partial differential equations as a series 
of ordinary differential equatiors in one of the independent 
variables is described. Application of a quasilinearization 
technique coupled with integration using Chebyshev 
polynomials to the incompressible momentum equation is 
outlined, as it eliminates the need for shooting to overcome 


the incomplete initial value set. 


An extension of the quasilinearization with Chebyshev 
integration technique to the compressible equations is 
proposed. BLi os control in view of the difference- 
differential scheme is investigated, and the difficulties to 
be encountered in applying the method are discussed. A 
comparison of the conventional and proposed approaches is 
made and the proposed method is shown to be potentially 


superior. 
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LABLE VORP RSYHBOLS 


SYMBCL DEFINITION 

Pp pressure 

p density 

KB viscosity 

Pr Prandtl number 

H enthalpy 

Ppp 
subscript e just external to the boundary layer 
subseri pty co remcte from the boundary layer 

Ab see [5] 

P (x/ueg) (due/dx) 

R (x/t) (dr/dx) 

N (P+1) /2+R 

By multistep integration constants 


Other symbols are defined in context. 
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The equations which describe the flow of a fluid must 
express the characteristics of this particular kand of 
physical problem [{10, 17, 19, 26]. The cohesive nature of 
the fluid is expressed in the equation of continuity; a 
volume element in the fluid must have a balanced transfer of 
mass through its walls if the molecular forces between its 
particles are to be obeyed. Forces acting on the fluid, 
whether external pressures or internal stresses, must find 
expressicn in an equation of motion developed from Newton's 
Second Law. Such an equation is usually given in terms of 
momentum as the velocity of the flow in different regions of 
the fluid is the quantity cf prime interest. If the fluid 
is incompressible, kinetic forces suffice to describe the 
system; se ag is compressible, thermodynamic forces 
resulting from volume and pressure changes must be taken 
into account. An equation based upon some convenient 
measure of energy must be developed from the fundamental 
thermodynamic relationships. This latter relationship is 
known as the energy equation. In addition, if other 
phenomena such as mass transfer through a boundary of the 
fluid are present, further equations must be formulated to 
buimdgcthemtsainto sacccunt. The equations to be studied in 
this paper arise in the study of the boundary layer around 


an axisymmetric body immersed in a compressible fluid at 


suue Sai2 = Io vans an? sulsonen a sencsiaye te 
to botd saAtwoheing eid? Yo soltebraraasado ate 
Ss ‘sduten svivedos dP sped el OF eof} wsldeta Lact 
@ eyiitarsion in ool Toupee Mat a beseatzgxe at Studs ca 
% te¥aner? baondiad ¢ sve) Yewe bint? od* ai taawets pavtov 


ati neswyod seoto! calodatfow odd SS elinw 222 dpvord? aso 


ed 
a 


ne 


vblat? oft go Babine aeszot, -sbeysdo S¢ of e378 aefobsaeg) 
fokt, Tale ,aseterze jkirpbulls to/ bebheng ivnistee tsdtou@ 
etaotvet? act? bogolevsh cuisom Io dedi supe be ar sabenoagxe 
to east ‘nd torte yllewad ed wottedps op foe = +Med baoped 
So @Qnotoet trst5t27b ai wold aft Io Ys toot ev aay es sptnowo 
Siuft ed? tI .teersgud oddug. le Pwiuidoy sit af bore 7) 
eit sditseck oo shbtiea #end0% obfeakn abitome<inel ale 
ge3G2 Siasaybouvads § ,sidiesezgaéos 2h #3 ws rmosnye 
naes> ad f2'n eoormeis sressty bus sapiov obits postings 
,Faeteguvio> ees aogn o¢ea4 a9rsetpa ‘1A n0o708 ‘etek 
‘tetnenshovt ait wort Beybleveh sd tenn (pteae Yo arueses, 
el qive@uolrntes wettsi eid? .aqadanmeysivt >tnany coms 
‘aerizo Ri wloieibhhs «st © .notrsups: eran. ar 48 aeowD 
aud 4c ¢anbnpod & apenas ed istanay? op tua “s 
ee Lclacusieeal as stoned snoksonge a sett 


hypersonic speeds. The boundary layer is that region in 
Which flcw perpendicular to the body walls is present. 
Although such equations are representative of a large class 
of prcblems they are apparently incapable of formal 
solution. Also, because of their nonlinearity, numerical 
techniques with well-defined systematic error control are 
avehicuLtgetot defines forththeerpreblen. However, their 
practical importance makes the development of methods of 
approximates solution vital; both the engineering design of 
high speed vehicles and also the extension of knowledge ina 
difficult scientific field are at stake. Challenges such as 
these have provided the motive for the development of 
numerical analysis thrcughout the history of that subject, 
and have laid the groundwork for most theoretical 
conSideration of the validity of solutions of mathematically 
complex problems. It is in view of the importance of the 
interaction between an experimental, practical situation and 
the stimulated theoretical investigations that this paper 


has been written. 


An effective technique for this particular problem was 
euclined tires by Clutter and Shathean) 196s ) [5] ) and» “other 
authors have made contributions since then [11, 12, 13, 14, 
20, 29]. The problem will be described and the various 


approaches to its soluticn discussed and contrasted. 


The standard description of the flow around a body 


immersed in a fluid consists of the equations of fluid 
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continuity, mcmentum conservation and energy conservation. 
Included in the equations are terms which describe the 
physical characteristics!) of the fluid: here, these include 
the compressibility of the fluid. Transverse curvature of 
the bedy is aiso included; such higher order terms are 
usually left out when Prandtl's boundary layer equations are 
used. A brief outline cf the derivation of these equations 
is given in Appendix A. The symbols used are defined in the 


Table of Symbols. 


A two-dimensional coordinate system is sufficient to 
describe the problem as it is axisymmetric [28]. Steady 
equilibrium ficw about a body of revolution is assumed. The 


coordinate system to be used is shown in Figure One. 
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One coordinate, x, is used to represent distances from 
the leading edge, or stagnation point at the nose, directed 
along the body surface; the cther, y, is used to represent 
distances from the body surface directed into the boundary 


layer perpendicular to x. 


In terms of this coordinate system, and with respect 
to the various fluid properties and thermodynamic variables 
described in Appendix A and the Table of Symbols, the 
equations are as follows. The dependent variables u and v 
represent the flow speeds in the x and y directions 


respectively. 


CONTINUITY: * (5. (pt) ie 35 (rpv) ]g= 0 
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The equations are awkward to integrate in this form, 
because the boundary layer thickness varies greatly with x, 
and numerical solutions will be seen to depend on values of 
the unkrowns and fluid properties at previous positions in 


the x direction. Moreover as stated the equations exhibit 
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singularity along the curve x=0 [5]. A technique which is 
often used ai Hydrodynamics is employed; a further 
transformation using potential functions combines the 
continuity and momentum equations while increasing the order 
of the problem [19]. A description of these transformations 


can be found in Appendix B. 


Lneetinal form of the equations ise 97 7 ee 
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Where X, A are independent variables ana, are 


dependent. pi’ are abbreviations EOE the parities: 
derivatives of the quantities with respect ton » Other 
quantities are physical parameters which are supplied in 
tabular form or calculated with respect to position. C, Pr 
and p depend on the velocity distribution (the solution ) as 
well as on position and must be included in the step-by-step 


solution of the problem. 


Bcundary conditions are: 
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Thus the problem has a free boundary with asymptotic 
approach to the outer boundary values of the unknowns as the 


flow merges with the uniform external flow. 


Two properties of the problem in this formulation 
should ke noted. First, the momentum equation involves only 
cne dependent variable, but the set of initial values is 
incomplete and the equation is nonlinear. Second, the 
energy equation involves both unknowns and is linear in p ° 
It is on these two features that the method of Clutter and 


Smith is structured. 


In setting up the discrete approximation of the 
equations the Hartree-Womersley technique [9] is used. That 
is, derivatives in x are expressed as finite differences but 
at each value of x the equations are treated as ordinary 
differential problems in q: The advantage of this treatment 
is that the extensive work dcne on the solution of ordinary 
differential equations can be employed, and as the equations 


are nonlinear, such an apprceach is doubly valuable. 
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The ‘following, is an outline ef Clutter dnd»Smith's 
technique: 


At a given x station: 


(1) Assume fluid properties Cc, Pr, f% from previous x 
Station (initially given at x=0), or from latest soluticns 


available for this value of x. 


“ 
(ii) Use a shooting technique to obtain 0 &Jana solve 


the momentum equation as an O.D.E. 


(a) Solve the energy equation as a linear problem 
in Vy) using the values of d, 9’, and gp" obtained from the 


solution of the momentum equation. 


(bo) Update fluid properties using latest values 


of all unknowns. 


(c) Repeat (a)-(b) until convergence is obtained 


on fluid properties and up or iterate as desired. 


(iii) Repeat (i)-(iii) until convergence is obtained 


on all quantities or iterate as desired. 


The procedure involves certain numerical components. 
They are identified below together with a statement of the 


chosen method cf attack cf Clutter and Smith. 


(a) Finite difference representation in x (third-order 


implicit differences). 
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(b) Integration in q (a series of Taylor's formulae of 
increasing step size to get started, then fourth-order 
Falkner-Adams predictor-corrector steps of constant step 


Size with one use of ccrrection). 


(c) Shooting tc get a complete initial value set for 
the momentum equation ( Trial and error until three "close" 


trajectcries are obtained). 


(d) Solution of the energy equation as a linear 
problem ( Superposition of solutions with integration as in 


(b)). 


(€) Choice of starting values of the fluid properties 
for this value of x ( the previous station's values are 


used). 


A detailed description of the technique is given in 


Chapter Ten. 


A general comment as to the nature of the Clutter and 
Smith method can be made. It is a direct approach which is 
built on the nature of the problem at hand and which makes 
use of the available information. Its value has been proved 
experimentally in a ‘Variety Of situations [2,57 11, 12, 
13,28,2¢]. ‘Theoretical consideration of the method was done 


after it was found to be effective. 


Since this basic structure waS proposed various 
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alternative treatments of components of the method have been 
suggested. in’ the folloWing chapters’ Clutter and Smith's 
methods will be outlined in some detail and contrasted with 


these alternatives in terms of numerical characteristics. 


The purpose of this exposition will be to illustrate 
fhesnuwesical) conplexity*tofea problemeroGelthis®? diftiiculty, 
The effectiveness of a technique developed under conditions 
of engineering necessity must not be denied; rather the 
desireakility of a priori knowledge of the behavior of the 
numerical components of a method will be emphasised. 
Finding an experimental method which gives experimentally 
verifiable results is an art; developing technigues within a 
context of knowledge of the homomorphisms between the 
mathematical spaces of exact and approximate solutions is a 
science and can therefore be used with confidence in a 
variety of situations and under known criteria for validity. 
The more complex a problem, the more necessary it is to have 
a guarantee of the behavior of numerical methods, and yet 
the harder it is to obtain one. A method which offers a 
degree of confidence which appears to be as high as may be 
obtained will be discussed. The application of remarks. on 
the numerical solution of the problem under discussion to a 
wider class of problems cf some similarity will be 


emphasised. 
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The Hartree-Womersley technique [9] consists of the 
replacement of stream-wise (x) derivatives by finite 
difference approximaticns tc produce ordinary differential 
equations inf). Derivatives are approximated at x=x(n) 
using quantities previously calculated upstream at x(n-1), 
x(n-2), etc., so the errer of approximation is cumulative. 
There is no difficulty in starting the solution as terms 
involving x derivatives disappear at x=0. At x=x(2) a two- 


point fcrmula must be used. 
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Two-point formula for function f: 
(o£) = fy 7 fw-2 x *N-1 
ax’ N (x (N) -x (N-1) ) 
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At x=x(3) a three point scheme may be initiated: 
Lagrangian three-point fcrmula: 
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Clutter and Smith decided to use the three-point form 
throughcut. The three-point form is contrasted with a four- 
point fcrmulae of the same type and with what is called the 
‘mean form’ using average values of derivatives over a range 
of points to estimate the derivative at the next point [28]. 
Tt was found that the point form is more desireable than the 
mean form as the latter process tends to diverge wildly as 
the step size h decreases. The four-point form, secondly, 
was found not to give significantly more accuracy for the 


additional computation and storage involved. 
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Consider the constant differences formula (1). Under 


this approximation the mcmentum equation becomes 


(MB ath), Py Myrna) 
(2 4B Baa) 


c, x(N) [(p'+1) : gn 


3 " é ‘ eee 
= an (Cyn) + T Go CyON - CAP By” + 20,41-—21 
+ Cy N(Bitn) Be Bk ee oD oe) me Oe OF 


Terms (10 Ow ’ One and Dna can be seen to involve 


the multiplicative factor x(n)/h. 


Clutter and Smith find experimentally [28] that if the 
ratio x(n)/h is bounded above by 25 the accuracy of the 
method can be sustained. However, a different approach to 
the problem of truncation error growth is taken by Jaffe and 
Ckamura [12]. They use one two-point and one three-point 
formula to start their three-point formula is a linear 
least squares approximation). Thereafter they proceed with 
a four-point quadratic least-squares finite difference 
scheme. Under their formulation the left-hand side of the 


momentum equation becomes 
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Py = Dece a: 2By x (N) ) where 
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hand side of the momentum equation involves the 


multiplicative factor x(n) instead of x(n)/h. 


Error analysis of a theoretical nature is impractical 
due tc the nonlinear form of the momentum equation. 
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CHAPTER THREE 
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In the Hartree-Womersley technique [9] integration 
with respect to is done along lines of constant x value. 
The x derivatives are approximated by finite differences. 
The equations are dealt with in a form in which the 
thickness of the boundary layer remains nearly constant with 
respect to the variation of x, that is, along the body. 
Under the transformation given in Chapter One the range of 
the integration variable n is [€,0) and the method of 
Clutter and Smith depends on choosing a finite final value 
TOL eS) eC tiatethe range: 1S approximated by [Oe%o) « Such 
choice is done by experimentation with the method and 
comparison of solutions generated with results obtained in 
wind-tunnel studies. Then the need for timing decisions to 
change step size and the associated use of interpolation on 
starting values of fluid properties are eliminated as the 
Tange of in 1s" constant. Clutter and Smith [5] decide 
therefore to apply a multistep extrapolation-interpolation 
(predictor-correctcr) scheme suitable for this type of 
problem. Allied to the use of this kind of integration is 


the establishment of starting values for the process. 


The integration problem can be restated as follows: 
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MOMENTUM 3 " p " 2 p 
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The procedure selected by Clutter and Smith to 
integrate these equations is outlined by Collatz [4]. It 


consists of the Falkner extrapolation formula followed by a 


Single use of the Adams interpolation formula. For reasons 
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of desired accuracy and speed discussed in [5], four point 


forms were used. They are: 
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of the Falkner-Adams scheme is to be preserved at the 
beginning. A detailed description of these starting 


formulae and the accuracy requirements can be found in [5]. 


The stability and convergence of the application of 
the predictor-corrector pair must be established. A 
discussicn of the derivaticn of the formulae can be found in 
{4}. An interpretation of the error growth is given by 


Collatz and the applicable theory restated here. 


Error bounds can be established for a repeated series 
of formulae if the prcblem being solved satisfies a 
Lipschitz condition on the convex domain of solution D. The 
nonlinearity of the problem and the nature of solutions for 
some boundary values will actually prevent the problem from 
Satisfying a Lipschitz condition in most cases; however the 
method cf error prediction used in this type of numerical 
solution will be outlined in order to illustrate how and why 


aQ.fficuttpes “occur. 
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solutions and approximaticns to those solutions, an idea can 


be obtained as to the rate cf error growth. 


Repeated uses of the corrector can be shown to 
converge under the Lipschitz condition and by suitably 
restricting the step size. Clutter and Smith point out that 
on e€xperimental investigation one use of the corrector leads 
to convergence within their tolerance at a reasonable 
stepsize. They choose therefore to ignore the problem of 
convergence. It remains tc investigate the stability of the 


process and the growth of truncation error. 
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different order derivatives). This system will provide a 
systematic error estimation process if the coefficient 
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because in the case that it is monotonic, the solution 


can be kEounded above. 


That is, if A represents the coefficient matrix and if 


an apprcximaticn to can be found such that 
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Then we can write 


(m) | ,n-m ™ ey is 
Ye] h b RPn-mjp ~ r+l a 
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n-m-1 * 25 

h (m+v) yn-m a (v) 

E = Y + 
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Such an expressicn is tedious to apply recursively. A 
more convenient guide tc the growth of error can be obtained 
by appiying a standard technique tc solve the difference 
equatiors for the error terms. That iS, a particular 


solution is added to a homogeneous solution of the form 
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The region where errcr growth is bounded is the region 

where these equations have a nontrivial solution. This 
region is defined where the coefficient matrix for the 

vanishes. ip 3 the equaticns under present study are 


concerned the system will take the form, where n=3, p=3, 


m=O, 1, 2 
p 3-m 3-9) 
Cc Ze = (h yee Coys (eee re: | Z 
z—-m Vv 
2 h = 
a vt “mev : 
v=0 


New as stated in terms of backward differences in 
their early paper [28] Clutter and Smith are using the Adams 
interpolator and the above theory is applicable. in, «their 
Fatereupaper -[5]}, atvmethod” closely gelinked soto) thecidans 
mterpolatorehut inet Gidentical to iteisestatedsinetterms T6f 
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the Lagrangian form. After the calculation of ¢ in terms 


c+\ 
cf Onc the lower order derivatives are calculated in 
terms. cf gp" instead of Ors, A study of the error growth 
in line with that above for a second-order problem could be 
done instead of for a third order problem consistently, but 


thé sLipschitz.condition would then have, to»stbessapplied sto 


(Ou which quantity is to be calculated at each step. 


The conditions for the determinant to vanish in the 


second-crder problem are 
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and the roots z can be found as functions of the Lipschitz 
constants ae K, £Or Q, Q’. The general solution of the 
difference equation for as will become dominated by the 
root largest in absolute value as r becomes large, and thus 
this roct provides a guide te the growth of error. Such an 
analysis using the formulae of Clutter and Smith takes no 
account of the bcundedness of tiber however, nor can it 
govern the error growth in Q at the first step in the 
interpolation. All that can be said is that assuming the 
existence of K, throughout the process the error in w) and 
(J will grow approximately as me where Z is the largest 


root of the akove process. 
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It should be neted that in regions where there are 
Singularities in the flow, that is where fw) and oy become 


unbounded, the error will have unbounded growth as expected. 


It has been found that the solutions can be extremely 
sensitive to initial conditions for some values of physical 
parameters in the equations. This is the problem of 
inherent instability of the differential equations. In such 
cases the approach of Clutter and Smith to the problen 
breaks down as the shocting method needs final values of gp" 
within some finite range to aliow interpolation. The 
Situaticn is aggravated if small increments in x are desired 
(see Chapter Two). The most effective approach to the 
problem of instability of the equations in such cases seems 


to be linearization of the equations; this approach is 


discussed in Chapter Six. 


The question as to the presence of partial instability 
in the finite difference equations is not so easy to answer. 
Clutter and Smith find that with too large a stepsize with 
some flows they are not able to obtain meaningful results. 
They attribute this to growth of truncation error. 
Certainly according to their results there does not seem to 


be a sudden threshhold in increasing step size where there 


is sudden growth of round-off error. 


There appear to be two investigations of alternative 


methods cf integration of the equations. Jaffe and Okamura 
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[12] suggest the use of Runge- Kutta integration for third- 
order differential equaticns. Jaffe and Thomas {14] use the 
radically different approach of Chebyshev integration on a 
linearized version of the problem. The latter technique is 


discussed as a separate topic in Chapter Six. 


Runge-Kutta integration is well-known. The procedure 
used by Jaffe and Okamura is taken from Collatz's book [4] 


and is outlined below. 
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The main advantage of Runge Kutta integration is that 
stepsize can be changed easily without interpolation on 
previous steps as it is a one-step method. According to 
Jaffe and Okamura this property is valuable as even under 
the transformed version of the equation the thickness of the 
boundary layer begins to change significantly in flows where 
integration proceeds far downstream. They can also 


economise on calculations by making use of the fact that as 
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increases away from the kody changes in the dependent 
variables become less rapid and therefore the step size can 
be increased. ite “their. applacetion of , the eiintegration 
formulae they include a means of testing these rates of 
change as integration proceeds at each x station and 


automatically adjusting the step size. 


There are two main disavantages of the Runge-Kutta 
type of process. In the first place, there are frequent 
function evaluations (in a fourth-order method, four to each 
step). Secondly, there is a lack of systematic truncation 


error estimaticn. 


Some indication of the size of the truncation error 
can be cbhtained by performing two integrations of different 
step size in parallel and using the fact that the method is 
of fourth “order. Over a large number of steps it may be 
assumed the rate of change of the error growth is dominated 
by jesse of the highest order derivative. THenG lsat 
integration is performed with step sizes h and 2h the 


follcwing estimation may be made. 


At each step, an error of order h* is made (agreement 
with aoZaylor®s series up to’ the “fourth order®*térm. is 
assured, where the expansion is around the previous point). 
At the point n=n,+2nhreached in 2n steps of size h the error 
Of truncation is A (anh? where A is some constant. Then 


at the same point Not anh reached in n steps with step size 
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2h the error of truncation is A (n) (2h) ‘ Then, if the 


process is being reasonably accurate, 
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where Pi da oe Y2LhyiL are the approximations obtained by the 


two integrations. Now 


(A) (2n) (h°) = A(2nh)h4 = atn-n,)ho 


4 
(n) (22n) (h?) = A(2nh)2°h” = A(n-n,) 2"h 


so 
A(n-n,)h* = @h,2n - g" 
a0) J EP Agee 
oP ae 
and 
= oan 
Bre cen — 2m n 7 Pn ,2n 
Thus the error at each value of r should be 


approximately one-fifteenth of the difference between 
integrations performed with step size h and step size 2h. A 
check on the reliability of the process can be maintained 


throughcut by maintaining two concurrent integrations. 


The use of Merson's variant of the Runge-Kutta process 
may be considered. The trouble is that when nonlinearities 


are present, the process of performing one extra calculation 
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in each step and using a comparison between it and the 
estimate of the point cf interest to estimate the error has 
been shcwn to consistently over-estimate the error. This 
drawback makes the method unsuitable for the problem at hand 
where changes in step size should not be necessary 
frequently and reaction to changing conditions should not be 


over-sersitive. 


The remarks that have been made about integration 
procedures in this chapter have centred on the treatment of 
the momentum equation. The same methods of integration have 
been used in most cases for the energy equation. In any 
case, since the energy equation is linear, the success of 
the process will be governed by the success with the 


nonlinear momentum equation [28, 5}. 
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CHAFTER FOUR 
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The solution of the nonlinear momentum equation is 
complicated by the fact that the initial value set is not 
complete. 0,9 is not specified. The technique Clutter and 
Suith [5] introduce to deal with this situation Consists of 
a process of shooting to obtain bracketting values of one of 
the final conditions within a certain range and then using 


interpolation to ccnverge on the correct trajectory. 
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Problems arise in the application of this method. In 
some cases, small changes in the starting values for trial 
trajectcries lead to large changes in the resultant curves 
{12,29], as can be seen in FigureM . Sometimes these 
changes are unmanageably large due to instability of the 
differential equations with certain values of the physical 
parameters. Even when the equations are not so unstable the 
difference between bracketting solutions may be too large 
for the Lagrange interpolation used to be meaningful. 
Ciurter and Smith [29] try “0 deal with this Slatter 
situaticn by stopping integration before the trajectories 
diverge too much, and repeating the process in the hope that 
the trajectory will gradually extend outwards to the full 
range cf rn as the correct solution is approached. Ina 
later modification the use of ‘mid-course corrections' is 
proposed as trajectories start to diverge. The problem is 
then treated as a new initial value one and interpolation 
done on the final trajectory segments of the resultant 


bracketting solutions. 


Jaffe and Okamura [12] attempt to improve the 
situaticn by introducing a form of quadratic convergence to 


the final trajectory. They introduce the function 
v 


(] v 2 2 tee 
at integraticn step n, where W, W, and W, are weights. Now 
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lies ne 
at the edge of the boundary layer o is explicitly zero and 
the cther quantities implicitly tend to zero, so W is made 


large in comparison tc Wo and W;. Then the minimum of the 


Sit we sade 
iaiis zat asutet pat? asim ont at eacaeee ite’ eases 6 


seetos Anarllaza: oF af aes (a spgel of pee t. zeit cet et: 
is ae 
escdd coufteec® . Deaneli cf nase sa M85 Sh 


ead5 302 Viet eIen, oy op 9@2bi yidesosaange ome a 2 


. - 133 

on 764 ag: > Vim 2i 215 oT ‘sy>ead teu Fad sos 78 ee 

i —, 

sib tente 7 ij J I , ne a) _ ie sb0s3pet an2 z¢ , 
oe a 

= ‘ - a - 

tc 4a9 ea [.? Ps a». y 77 } 15 Fé G GOA - : in £2 


asitorversid 4)'7. 210% ageTratbsTa. nrarauvere va 


ads ays (A j 9 > ply! : nus Dui out aprevi 


>. > 
I0n* e=> 7 = ~ =% ’ ‘ vg _) a 4 I jeg ad9 
é , 
Saotvures wmerts 
Fs 
. Ae : = 
eat. efome: of vrypsers (07) Smyesdo tgs BR 


OF S0uuPsNYReD S2767hsUG Io: 2152 eniowigsenh qa. wi 


aétbisanl siy =spboysn: ven 
, - j 


32 


function £0) is equal to zere, ‘and the? location) "of! *this 
Pee ss i / i Md 
Minimum will yield the desired values of Oe ana D ° 
1 
Then p'K)9 can be improved to yield these values by a process 
of quadratic interpolation. The authors find that this 
a . / 
process 1S rapidly ccnvergent for problems where “) is such 


that it is monotonically increasing to zero. 


Any interpolation method is dependent on finding 
bracketting values of the desired soluticn within a 
reasonatle range; the variations proposed can accelerate 
this search but they cannot eliminate the cases where 
inherent instability is present. It is to this end that 
Jaffe and Thomas [14] developed the idea of linearizing the 
equation. The sclution of the problem in terms of the 
difference between succeeding trajectories in a quasi- 
linearization iteration means that smaller variations in the 
numbers ccncerned may be encountered in situations where 
instability is not extreme. However, once the problem can 
be treated as a linear one the need for shooting can be 
eliminated entirely by treating both endpoints as fixed and 
solving the trajectory as a system of linear algebraic 
equations. The one conditicn that has to be met is that the 
initial trajectory be ‘sufficiently close' to the final one 
to ensure convergence of the quasilinearization process. 
Jaffe and Thomas appear to find this condition easier to 
meet than it is to overcome the problems involved with the 


more conventional approach of shooting. The 
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CHAPTER FIVE 


The energy equation can be treated as linear problem 
« ! iu 
0 \) alcne once od, 0 and Q are considered supplied by 
the current solution of the momentum equation [28]. The 
integration formulae are discussed in Chapter Three. The 


boundary conditions are 


n = 0: either y(0,x) = g(x) - 1 
OL ew Oye) “Seq tise) 


ie ere yp > 0 or w' + 0 


Linear problems of this nature can be solved by 
constructing a linear combination of the solutions of two 
other fproblems. Clutter and Smith [28] and later authors 


(11,12,13] use this technique. 


/ 

Suppose UR is given. Integrate with arbitrary VA = Uy 
to obtain UU ied y' is greater than zero at the edge of 
the boundary layer, subtract a quantity from the arbitrary 
starting condition §to obtain Wa and integrate to obtain 


ye Then calcuiete AB such that 


Ay, + BV, = Av, + BY, = 0 


at the edge of the boundary layer, and 
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Then since 


= = = O 

bean ty.2 = hy 299 
then §=AU, + Bb, 
and j. = Avy fe Bys 


will satisfy all boundary conditons. 
A parallel procedure is followed if Whdis given. 


Once both Vy and w) are known, fluid properties along 
the constant x line can be updated. Iteration on the energy 
equation can be performed improving these fluid properties 
until desired convergence criteria are fulfilled. Then the 
newest fluid property values can be fed to the momentum 
equation and the entire doukle loop perfcrmed again until a 


final set of convergence criteria is met for the x station. 


Various authors [11,13] have inserted other linear 
equations describing phenomena in special flows (for 
example, gas injection into the boundary layer), either in 
place of or coupled with the energy equation. The treatment 
ef such equations is identical to that of the energy 


equation. 


Any technique for solving "linear" ordinary 
differential equations could have been used in the solution 
cf the energy equation. These techniques do not affect the 


fundamental nature of the Clutter and Smith [28,5] method 


Set panic bas aeysuto ee? a0 autem Luba 


was 108 % 


2a0tibdos yasbanod Wee 4 
eh 
: 7 ’ 

. pi : 
sovin arhayy ++ fevol fet 2: s2ubss0rq letlerta ry yy 
| 7 

pe 
‘+ 


gftele sets seesrq Si61S \ueony ers O tae Y. dod sand 
yerecs sils 29 a@apaies)). -haiedgy Sa” dso ecki = oa 


dekeveqerq wint? sapd? ypelypr jaf feertetzeg ec 780 


SA¢ ned? .bef{[ltivs are siasziip epaspis¥ aps bo31F256 
qigeedd 237 5S Sot od. ase, 25dfayoyseagortd feat 
@4igay ateee feel 5T34¢ goo selineh exiete als bas 


swO2ente x ety 10% ged 20622597245 Bouspisvass Fo fae 


d6eabi -ae4d2e batsceni aved fet fh] esattrs enotrey om! ; 
$53) aWolt (sinsn= vi, shetpesdg ofideyneen cook tape” 
Mi gdetet> .(35yal yiahnved aft ores aetrostor aap signee 
feeesesct =dt «.tattsope vypvene ado driv Selgten go Ges 
Wesene edt Yo jsu¢ of sortaehl st =rpkrseupe Aaya 

«hods, 


yasatbae Waaenii" . ¢0iLvice 293 - Syeintoes 
ant rulos aft at Seek ASS 2454 fLluoD sin dasppe 
RNY #9ATE For Oh Soup tuiioss » 2/1 ha a 


v8 


it : 


a} 


36 


for solving the coupled equations, nor can a different 
choice cf integration method for the energy equation avoid 
the protlems that arise due to the nature of the momentum 
equation and its boundary values. A method of solution for 
both equations using a series expansion will be discussed in 
Chapter Seven, and this technique will be found to take full 


advantage of a linearized fcerm of both equations. 
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The use of quasilinearization in the approximate 
solution of differential systems originated in the theory of 
dytamic programming [1, 16, 25). The formulations most 
familiar to the numerical analyst are those of the Newton- 
Raphson and Kantorovich approximation techniques. The aim 
of the quasilinearization method is to set up an iterative 
process on a linear problem which approximates to the 
original nonlinear problem in the sense that the solution of 
the linear problem is a satisfactory approximation to the 
solution of the original problem. Then at each step in the 
iteraticn procedures suitable for linear problems may be 
used. The question of the validity of the resultant 


convergent solution must be investigated. 


Any high-order differential equation can be considered 
as a system of first-order differential equations. Consider 


therefore the system 
U = F(U,t) 


where U is the "state vector’ or set of dependent variables 
and in this case t is the only independent variable. Let 
the boundary ccnditions be 
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= 
Let U(1) be an initial guess to the solution of the 


= S| 
system, Then F( WU , T) can be expanded around U(1) ina 


Taylor series truncated after the linear term: 


F(a,t) = F(G,,t) + 9, (F) oly. (1) 


where 
a OF. 
wee . (sa) == 
2: U=U, 
where the indices I and j mean the ith and jth elements of 


the vectors F and U respectively. 


New it is desired to obtain a second approximate 


—_ a 

solution, U(2), which is closer to the real solution U than 
— — a 

was U(1). This can be achieved by setting U equal to U(2) 


in (1) and forcing 


then 
U. = F(U, ,t) + J, (F) (U,-U,) 


an iterative process can then be set up of the form 


Since at each stage in the iteration the solution 
depends on a two-pcint boundary problem, it must be shown 


that this sequence exists; if it does exist, it 1s necessary 
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39 
to examine its convergence and determine whether its limit 


is indeed the soluticn of the original problen. 


The linear iterative scheme may be written 


This is a linear first-order ordinary differential 


equation in the unknown U(n) which has the general formal 


solution 


where c(n) is chosen so that the boundary conditions are 


satisfied. 


max max(|1F]|, |lgi{) =m 


aes 


Diet 


and min iene 1 


then ies 


ly 
nee < exp (mT) [45 (m + m|{U,_411) exp (-m,T) dt + CJ] 


A sufficient condition for a convergent sequence to be 


obtained from this relation is 
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T max 


—\ 
HO. 1 < exp (mT) (2m. J, O<t, ST (exp (-m,t,)) dt +c.) 


exp (mT) (2mT + Cy) 
then if 


= 
1A 


ett om + Cy) 
for 


—-mT 
m. = 1% <e < 2mT + Cy 


thus, sufficient conditions for a convergent sequence to 


exist are 
1) m to exist and be bcunded 


2) U (1) to be chosen so that 
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emt 2mT + Cn 


The sequence can be shown to be monotonic as well as 


convergent. 


The preblem can usefully be rephrased in terms of the 
difference between succeeding approximations in the 
quasilinearization algorithm: 
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with boundary conditions 


and the above discussicn may be fcllowsd fer the unknown 
— 
aC). The latter treatment is valuable in investigating 
the rate of convergence cf the sequence and in determining 


Ene Limit point. Thus 
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= 
(exp 0 Jydt,) dt] 
where c(n)=0 as the boundary conditions are of the form 


e = 0 at t=O 7 ta 


The mean-value theorem applied in vector form yields 
=) eat — 1 
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where H (n) represents the Hessian tensor of second 


= 
@erivatives applied at a point in the anterval “[U{n=1), 
U(n) J (element-by-element). Now 
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Thus the rate of ccnvergence of the sequence of 


approximations is quadratic. 


By a process of inductive reasoning it can be seen 


that 
2 N 
= KT 2-1 = 2N (2) 
llega tI acre exp (mT) ] (max|le,11) oe 
where |. es an 
Soy a 
Cenvergence therefore depends on the quantity 
2 
KY max mane 
gp oP ME) jeepers) Goh 
This quantity can be made less than 1 if either T is 
. Ba 
small encugh or sere lleo1| 
is small enough. Thus a proper choice of the initial 


approximation can be sufficient to ensure a quadratic rate 


of convergence. 


If the ccenvergence condition above is met, then the 
Sequence of approximaticns 2 (n) has. Limit eporoteG. tit 
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The quasilinearization process is particularly 
appropriate when used in conjunction with the Hartree- 
Womersley technique [9], because at each x-value except the 
fanst an initial approximation is available from the 
preceeding station. Not only that, but the initial solution 
is known to be within a certain tolerance of the actual 


solution dependent on the x derivatives. For example, 


3 2 
(P) cex(m) = )yex(m-1) * (32) (x - xy) + Olax)* (3) 


m-1 


and similarly for the other unknowns. 


What this means in practice is that a bound can be put 


on the quantity 


a ay = —_ | _a 
Ure Oy ile ed ater eel, ees by monotonicity 


It would be possikle, indeed, to ensure convergence by 
keeping running estimates of the relevant quantities in 
equation (2) and calculating as many terms as required in 
equation (3) +. Coincident with a decision to keep such a 
running check on convergence would be a decision as to the 
appropriate kind of norm to use. A situation which would 
make such a precess valuable would be, for example, the 
approach to a point of separation in the boundary layer 
equations. The x-array could be manipulated as well as all 
the other quantities involved in the ccnvergence criteria in 
order to guarantee convergence and the validity of results. 


Problems with the use of smaller increments of x as 
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discussed in Chapter 2 are avoided when Chebyshev 
integration is used, showing again the power of the combined 


techniques. 
The running calculations involved are as follows: 


(a) suppose H and J are known on line oo. Estimate 
them on line m+1 by element-wise Taylor expansion. (since 
the elements of these matrices will consist of multiples and 
combinations of the coefficients and unknowns U this is 


possible). 


(b) 

Ensure 
max 2 aly a ol 2 ene 
O<t<T llegll < 110 U,1l = =a exp (-mT) 


by © "choice ‘of 0 (1) threugh equation (3) and/or choice of 


x (ajimx (a 1). 


Once a quasilinearized form of the equations has been 
obtained the choice of a method of integration has to be 
reconsidered. The incomplete initial value set remains, and 
unless both boundaries are treated simultaneously a form of 
shooting must still be used. The use of the difference 
between trajectories as the dependent variable will reduce 
this problem as variation between different solutions for 
different initial values will be less, but the effect of 


this reduction is hatd to predict. 


A method of solution of a linearized form of the 


incompressible momentum equation has been described by 
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Cebeci and Smith[2]; ite uses /a"finite difference scheme 
designed to take advantage of the nature ‘apa the 
incompressible problem and therefore has no application to 
the general problem. Other experimenters such as Radbill 
and McCue [25] describe the construction of solutions to 
this type of problem using any integration method for 
ordinary differential equations and superimposing 
homogeneous and particular solutions. The latter approach 
loses accuracy when the superimposed solutions are greatly 
different for small changes in the initial conditions; it 
has been demonstrated that this is precisely the difficulty 
with the present problem. It is imperative, therefore, that 
a form of solution which treats both boundaries 
Simultaneously be selected. Such a requirement suggests a 
means of transforming the differential problem into a system 


of linear algebraic equations by means of a series solution. 


The application of quasilinearization to the boundary 
layer equations combined with Chebyshev integration will be 
discussed in the following chapters. It should be added 
that in general nonlinear bcundary conditons can be treated 
in exactly the same way as the nonlinear differential 
problem. Convergence to the actual boundary conditions 


follows in a similar manner. 
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At each step in the iteration of a gquasilinearization 
process integration “of ‘a’ vlinearc hordinary® differential 
equation must be done. Although an analytic solution was 
presented in the previcus chapter in crder that convergence 
analysis @enighty beotdone, “fhiss *soluti6n Taisvenot useful 


computationally. 


In the problem under discussion one of the major 
difficulties has been the incomplete initial value set. 
Shooting is less risky using a quasilinearized form as 
differences between neighbouring trajectories are not so 
large when @ is used instead of v. It is more sensible, 
however, to take advantage of the linear form of the problem 
to the extent of using an integration technique which treats 
btoth boundari2s at once. Superposition of solutions as 
discussed in Chapter Five might be used, but once again the 
accuracy of the multipliers would be called into question if 
a degree of inherent instability were present. Jaffe and 
Thomas [14] have suggested a method for the incompressible 
momentum equation which converts the problem into a _ system 
ope innear algebraic equations, thus treating all boundary 
conditicns simultaneously. Other interesting features of 
the method will be discussed. Their choice of integration 
technique involves expressing the dependent variables in 


terms of Chebyshev polyncmials. In order to discuss this 
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choice cf procedure a brief outline of approximation theory 


will be given. 


Numerical analysis is a discipline which attempts to 
describe mathematical functions, which are defined either 
explicitly or implicitly, as a simple computational problem 
for various values of their independent variables. Such a 
description may consist of tables of values or, more 
economically, of an approximation function which is easily 
evaluated for given arguments and which has known deviation 


over the range of interest. 


A differential equation is an implicit means of 
defining a function or functions. Conventional step-by-step 
integration builds up a table of values which describes the 
function at discrete values cf the independent variable(s). 
Error in the table is cumulative, and the value of a 
Lunct#@on ateansarbitraty pointeanot: tin stheretables smusty pbe 
obtained by interpolaticn, wherein further error is 
introduced. These conditions are acceptable if enough is 
known akout the error. A different approach is to assume an 
appropriate approximate form for the solution in terms of an 
easily computed function of known deviational behavior and 
transform the differential problem into one of determining 


the parameters which will ensure a reasonable approximation. 


The easiest approximational form in terms of 


computation is the polynomial. The most compact ina given 
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Space is the orthogonal set of polynomials spanning that 
Space. Then, if a theory of' measure for the difference 
Fetween functions in the space and members of the basis or 
linear ccmbinations of members of the basis is defined, any 
function in the space may be represented non-redundantly in 
terms Of the basis «with knewn deviations Given measure 
definiticn, the linear combination of basis members 
Yepresenting a function may be manipulated in order to 
Minimize the difference in some sense. That is, if f is 


represented by 


pete S EAT) oh B50) tea eee 


the coefficients a(i) will depend on the definition of f and 
on some error criterion. The choice of norm will result in 


a certain set cf orthogonal functions being used. 


The Chebyshev polynomials are useful in numerical 
analysis because of their computational form, the error 
criteria they may be made to obey, and their applicability 
tomaewide ClaSs!or finctisnepielanCzosep J oqei Se ethesia cher toe 
modern interest in them. Since his re-discovery of their 
importance many researchers have investigated their 
application to approximation theory. A comprehensive book 
on their use is. by Fex and Parker [8]. In this book the 
means of extending known results to differential problems is 
discussed. The relevant theory and properties of Chebyshev 


polynomials is outlined briefly below. 
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Suppose that the Chebyshev polynomials are defined by 


T(x) = COS .((n=1), 6), ..cos_ Oi=]=.s, -c=1exs1 
or equivalently 


Ty = be T, = xX 
and the recursion relation 


Tey (*) ue De T(x) - 7T (x) 


bated b 

These polynomials obey the minimax criterion, that an 
approximation Pax) of f(x) , where e,(x)=£W)-p,%, should be 
such that max, | endl is a minimum for all possible 
representaticns in terms of combinations of polynomials. 


The Chebyshev polynomials may also be made to obey a least- 


squares criterion 


n+. 3 
a w(x, ) e7 (x3) = min, (discrete f) 
or 
: 2 
ey w (x) e, (x) dx = min, (continuous f) 


where the weighting functions are 4 
w(x) = (1 - x) 


Ancther property (ope the Chebyshev polynomial 
approximation is that, of all polynomial expansions 
expressed in terms of ultraspherical polynomials, it has the 


fastest and most universal ccnvergence rate. 


Evaluation of the coefficients for a function f 


involves expressions like 
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which, as given here, satisfy the least-squares error 
critericn. Where f is implicitly defined these expressions 
cannot be used. The coefficients must be calculated using 
other information. The differential system and knowledge of 
the differentiated form of a Chebyshev expansion may be used 
to set up a system of linear algebraic equations which may 


be solved for the coefficients. 


In order that the differential equation may be 
represented as a system of algebraic equations all 
quantities that appear in it must be represented as 
Chebyshev series, and means cf combining the expansions in 
products or quotients or in various other ways developed. 
This may not be a trivial task. The following chapter will 
illustrate the process as developed by Jaffe and Thomas 
[14], and the extension of the method to the more complex 
compressible equations will be discussed following that. A 


useful formula is that for the product of Chebyshev series: 


b n 
a n 1 
- go> pot) be 
eset to) a, T; Batiat ‘agisondh 
1i=2 
: Ee 44-9) Ob (i-3) 0] 
T.,T. = cos(i-1l) 0 cos(j-1) 0= 5 [cos (it+tj-2) G@cos (1i-Jj 
J 
af 
eet iyi, 
a,b n n a.b, n a, 
1 1 Lot 
£f.g = ( L + ae = a,b,) T, + ay eee ay Be geet 


pdiew faba tysieds sit vt nga en . 

to apietveny dns wataye Epidifsssat ib aes #7 eae 

fd20 od yse nclsseuts vadayaan? 3 ie atot Sereda aus odd, 

q2e. dnidw 2ucisevps sterdse.s in ees te esters oe tee 0! 
ehtsEoR3eso ent ae) 


sd yee woltevp> Lstrisre33fh. edt fede gebio 
lis ‘enotrsups 3ietdeole to gste¥e 8 28 
ae $52ne23itdor st seen, 3% gl qnoaiy's sade 
d* gacvkansaxe dt poinidgso to eddea Sos. ,jaeisee 
sPoqolsy¥a6 2yar seitto sugtisy “tito séeerdoae me! | 
tliv wrqgado poivoliol ai?  .s#e87 {fsiviat & st tow tne 
2s0ent Bus 232s "d Psqolevsb 28) 23530%9 eff ote: exh” 
xoiqué> o204 sit of Borten Sdt $6 aotensduo add bas (Pt) 


7 . - 


A ,28d¢ eaivolfo®? beeenoznh of IfLiw agcoftsepe sidiags gaol 


a 
— 


1. : 
tsealts? vederdar> 6 toubo7G 413 263 eadt et oe . bs) ix 


f 5 <2 ,6 ee 
+ :* ° aT - 


; - 
U6 (h-t)2o0% (S+E¥t}eee] 2 + (i-t}eos 8 (L-Red = 


er 7 ta é- 


1 


Various techniques for solving linear systems may be 
used tc obtain the Chebyshev coefficients. It should be 
borne ir mind that at least part of the system, that dealing 
with the differentiation cf a Chebyshev expansion, is upper 


triangular: 
/ Za 
ee SPs sor. (De oe ee) 


/ 


On = One 

The solution of a finite set of equations for a 
truncated Chebyshev expansicn involves a perturbed system 
and does not yield the leading terms of the infinite 
Chebyshev expansion. A methcd of estimating and partially 
eliminating the perturbation developed by Lanczos [18] and 
discussed by Fox and Parker [8] may be used but it means 
more calculation and in an iterative scheme such as 
gquasilinearization may not be valuable. The alternative is 
to calculate succeedingly larger sets of equations and 
examine the sequence of coefficients for evidence of 
convergence of the expansion. The latter approach was used 


by Jaffe and Thomas . 


A brief summary of the advantages of a coupling of the 
quasilinearization process with Chebyshev integration for 


the boundary layer equations will be given. 


A) Quasilinearization allows the problem to be treated 


as a linear cneée. 
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B) Quasilinearizaticn used in conjunction with the 
Hartree-Womersley technique [9] can ensure convergence to a 


solution where conditicns on the various derivatives are 


Sapeeried. 


C) The Chebyshev approach allows treatment of all 


boundary conditicns at cnce, thus eliminating shooting. 


D) The representation of solutions as Chebyshev 
expansicns is compact compared to tabular values and rapidly 


convergent ccmpared to cther polynomial expansions. 


E) Derivatives or integrals of given expansions are 


easily calculable. 


F) Any evaluation of a solution for a given value of 


the independent variable is possible without interpolation. 


G) Error bounds cn solutions are calculable. 
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In 1970 Jaffe and Thcmas [14] published a paper which 
solved the incompressible momentum equation using 
quasilinearization and Chebyshev integration. Their work 
serves as a practical illustration of the various procedures 
involved. If the equation is written as 
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CF Garay) 7 COQ may 9) Ogg -X yg?) 


then the quasilinearized form is 
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aE : eT x het (U5 West 
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Ccnvergence on the iterative scheme may be expected 
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subject to conditions outlined in Chapter Six. 


New suppose that the solutions are of the form 
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introducing (4m-6) unknowns. 
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yielding 3m-6 equations. 


Matching the coefficients for the polynomial T(j), the 
resultant system of linear algebraic equations corresponding 


to the differential system is 
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Soil ®) 
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yielding m-3 equations. 


where cheb, (F,_,) cheb, (F(U,_4)) 


cheb, (F(U;_5 + €3_4)) 
The number of equations is 3m-6+m-3=4m-9. The boundary 
conditicns give the remaining three equations needed for a 


completely determined solution and they take the form 


m 
e(-1l) = = - a, + a, - ...(-1) ae 0 
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e'(1) = 4 + ay + a, Sa eG an-l = 0 
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Scluticn of the system can proceed for each cycle of 
the quaesilinearization precedure given an initial solution. 


The follcwing initial solution was found to work: 
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5 
£o)=qen (2) [it wh ae) 
results of Jaffe and Thomas have been duplicated 
erder to correct their mistakes and as a basis for 
of the method to the compressible equation as the 
of solution is the same. It should be noted that 
Thomas obtain results which are comparable in 


and efficiency to both finite difference results 


and results cbtained using Clutter and Smith's [5] methods. 


The advantages of storing solutions in polynomial form have 


already been discussed. The following chapter will be 


devoted 


extending 


to a discussion of the problems involved in 


the technique to the compressible equations. 
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CHAPTER NINE 

The system consisting of the compressible equations, 
momentum and energy, may be quasilinearized in the same 
manner as the incompressible momentum equation. Indeed this 


has been done, as can be seen in references [1,25]. The 


equations may be written as 


Ng "' B" ,O' ,O,1', 7,4) = 


Pe 2 Pe 
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and then with the state vector 
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the quasilinearized form is 
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Ncw since the terms C, pep, T dre not Constants oF 
functions of position but depend on the solutions or state 
vector, the Jacobian matrix must take this dependence into 


account. For instance, 
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Thus either the exact form of the dependence of these 
guantities must be known or some other means of updating the 


Jacobian at each step must be found. An alternative 
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approach is to treat these quantities as being constant at 
each step and updating them after a new solution is found 
Without including the chain differentiation in the Jacobian 
matrix. The effect of such a measure would be to slow down 
convergence of the quaSsilinearization algorithm an 


unpredictable amcunt. 


A far greater difficulty in extending the Jaffe and 
Thomas [14] method to the compressible equations is that of 
expressing the fluid properties as Chekyshev series. As far 
as is known to the author this has not been attempted, 
although the desireability of doing so is great, in order to 


eliminate shooting. Problems arise, for example, in 
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can be expressed in terms of the Chebyshev expansions for Y 
and Os but an expression for the fractional powers of a 
Chebyshev series is not known. The above form of dependence 
of the fluid properties on the solutions was deduced from 


tables cf experimental figures by Clutter and Smith [5]. 


One approach to the treatment of Ce eok P., etc. woulda 
help in the calculaticn of approximations of the chain 
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expansien for ali desired quantities in the most 
straightforward manner at the cost of spoiling the elegance 
and suitability of the method and adding extra computation. 
However the ability to eliminate shooting is the paramount 


consideration. 


Suppose a new value of x is to be dealt with. Initial 
approximations to the state vector and fluid properties are 
available from the preceeding x station and they are assumed 
to be in Chebyshev form. Then one iteration of the 
quasilinearization scheme with the Jacobian matrix assuming 
constant fluid properties will yield a new state vector in 
Chebyshev form. Now new fluid properties in Chebyshev form 
must be obtained. The only possible way to do this seems to 
be to calculate the actual values of the relevant solutions 
at various values of ro» use Cohen's formulae to get the 
fluid properties at those points, and then fit a Chebyshev 


interpolation polynomial to then. 


The most accurate Chebyshev interpolating polynomial 
of degree n joining the points x(0), x(1), «++, x(n) is that 
calculated for x at the zeroes of the (nt1)th-degree 


Chebyshev polynomial T(n+2): The collocation points are 
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where 
a i fl 
app Rencren! #6174) 
1=0 
and 
ome 
oar = n+L ae f(x, ) Ty (x;) a= Fb Oy Flags” 


Uherverror- is then of order Tint2) % 


One fact serves to compensate partly for the labour of 
additional computation at ¢ach quasilinearization step. 
Cnce the fluid properties have been updated according to the 
new state vector, two succeeding expansions of the fluid 
proerties are known. Then an approximation of the form 
3B") Goo ic = C= Gy A Oy RG 
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accelerated form of the gquasilinearization scheme can be 
used. The above approximation can be formed, as can 
expansicns for quantities like T, 
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by use cf a quotient tormula fcr two Chebyshev series. 
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a Chebyshev expansion for 1h, in terms of the a(i) is 
desired. 
New 
1 
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the product formula for Chetyshev series yields 
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quaSilinearization steps. 


Under the suggested treatment of the fluid properties 
Jaffe and Thomas's [14] means of solving the incompressible 
equation can be extended to the compressible equations. The 
integration problen, with its associated error 
considerations, becomes a problem of solving a system of 
linear algebraic equations with associated questions about 
the well-conditionedness cf the coefficient matrices. The 
form of the matrix is given below, and as the elements vary 
at each value of x and each quasilinearization step the 
condition of the matrix has not been investigated a priori. 


If it becomes ill-conditioned the method will fail. 


The next chapter will compare the proposed extension 
of the Jaffe-Thomas technique [14], quasilinearization with 
Chebyshev integration, to the Clutter-Smith method [5] of 
solution. The design of implementations of the methods will 


be outlined. 
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CHAPTER TEN 


Flow diagrams for the two techniques are giver below. 
A list of the routines involved and their functions follows. 
The basic Clutter and Smith method [5] and Jaffe and Thomas 
routine [14] were implemented in order to compare 
programming complexity. These implementations were done 
using Fortran IV. Extensicn of the proposed method to the 
compressible equations is represented in the form of flow 
diagrams and descriptions of additional routines required. 
The three methods will be referred to by CS, JT and EJT 
respectively in the following comparison. Conclusions are 
to be drawn with reference to a class of equations bearing 
the characteristics of non-linearity and an incomplete 


initial value set. 


CS Routine Structure 
MAIN: Sets up arrays 
Sets N=C, M=0, E=0 
Sets initial conditions as necessary 
(a) Branch to treat x (i) 
increment I 
test for end of program: i=N 
if not, go to (Ca) 
Treatment Of X(i): 
MOM: Set anitial bp" for given trial 


integration 


a yy «2 Oy al 


ae 
Zaugragoar 38 wozies 
- . if: 
woled  asvip ote 2eppfadpet ovt od9 103 eusipstt vol 


‘ou a 
syaliot encrroaut aleds Bas bsvfovat zon sun adt 2: ell 


agmGdT Sire oF FAN Ons [OC ) bottom dslac ons zers0ls stead | sat 
S 

: ‘ = 7 aoe 

of 126706 rl Seatnosslaat sien [\0F] eat: yor 


a 


sich oiew anottstases([qeal szeqT “Yrs igeoo ecieeetQeag 
, 


Ay 63 often Ssagqozq Saf to noiensteS «VF neEdes t pater 
: ; one - 
vol® In e401 adt ai“Hstnseatess . ei <caoitasps oftteas 1#O% 


‘ * . = - 4 
batingst 2zsaidgbor Msnerstibs t6 gaotfaelaoret bas susxpadd 
a 7 - 


_ 
« 


"LS bes Th 22 vw of bovaete: od 4 


- 


a eholtse 67 ane dT 


“h ‘spoztauf9a63 ~Mectecgo> privollo? siz wr ylewe 7292 
} y . 


> B OF SodScetaa dzitv aveth Sd Gt 
7 9 
d 


stefqwosni is Bas viltsaail-to, 2 ‘eviteitesSareds sae 


vu 
+ 
i 
ad 
A 
i 
‘ 
tomy 
’ 


SYeL%s an eate2 
O=G .0=H , = etse 
YESEEeDGu es Sno liibno> fet ini aioe ; 


(£)x dseTt of fonstH (ap. 


Liaameazooe 
a=i hae tec thle toa 5) 
=f $@62~01d I6 Oas 192 +293 

- F 7 e . 


7 
B) OF op ,fou 3, 


> a 
i _ 
“ft { ni tted 2 


66 


Generate trial trajectory (IMOM) 
Use shocting logic with several trajectories 


(INT) incrementing M 


ENS Implement superposition logic 
Generate necessary trajectories (IEN) 

SErTrs Generate new fluid properties using 
given solutions and incrementing E. 

IMOM: Integration logic for D 

INT: Interpolation for shooting logic 

IEN: Integration logic for Y 
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COVAL: 


CHES: 


1S ONES 


FN: 


COBPE: 


sets up arrays 

Sets up initial approximations 

Sets parameters for the particular 
problem (e.g. Falkner-Skan) 

Calls COVAL to set up Jacobian 

Matrix for initial approximation 

Calls CLIN to implement method for each 
x value 

Evaluates Jacobian elements for 

latest approximate solution 

Sets up convergence test for line x(n) 
Replaces old solution by new solution 
Calls CHEB 

Calis ESOL 

Calls FN 

Uses Gaussian Elimination (for example) 


to solve system outlined in COEFF 
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Obtains expansions for solutions from those 


for the error terms 
Sets up the system of linear algebraic 
equations for the quasilinearized 


version of the problem. 


Obtain Chebyshev expansions 


for fluid properties using evaluation 
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NTS Interpolate with Chebyshev series 
Changes: 

MAIN: Different initial solutions 
COVAL: Different Jacobian elements 

QLIN: Extended ccnvergence tests 

COEFF: Different equations to represent 


quasilinear form of the compressible 
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Ccmparison of the methods is to be made on the grounds 
of thecretical considerations, programming considerations, 
computing time involved, and storage involved. 


1) Theoretical Considerations. 


SS cee Se ee ee See ae ee eS ae oe a ee 


(a) Both methods use the Hartree-Womersley technique 
to transform the problem into a series of ordinary 
differential equations. In the CS method, integration using 
predictor-corrector formulae is then used. First 
approximations to fluid properties are taken from the 


preceeding x value. In the 


ty 


JT, a quasilinearization 
process with Chebyshev integration is then used, with first 
approximations te all quantities given by the previous value 


Cleex 


(b) Integration. cs uses multiple individual 
integrations with different initial values and fluid 
properties. Each integration uses complex starting formulae 
and involves truncation error which cannot readily be 
estimated . Different multistep formulae are needed for 
each dependent variakle p* p', 0, oe / é A double 
iterative loop is needed to produce solutions for a given 
value of xX. EJT uses a Single iterative loop in its 
quasilinearization scheme , each loop involving solution of 
a system of linear algebraic equations, evaluation of the 
Jacobian elements, and peripheral treatment of fluid 


properties. 
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(c) Shooting . CS uses a shooting method to complete 
the initial value set. This process introduces many 
additional ccmplete integrations along x(“) and in some cases 
does nct werk due to the inherent nature of the problen. 
Suggested variants may improve the solution in borderline 
cases but de not eliminate the problem. EJT eliminates 
shooting entirely by transforming the problem into a linear 
one and using a method of integration which can treat both 


ends of the integration range at the same time. 


(d) Soluticn of the energy equation. CS uses a second 
and separate iterative lcop to solve the energy equation. 
In the EJT , the two equations are solved simultaneously in 


one iterative loop. 


(¢€) Fluid properties: the fluid properties can be 
treated more easily in their actual form in CS , but give 
rise tc more iteraticns. The interdependence of the 
solutions and the fluid properties is represented more 
accurately in the quasilinearization method, however, thus 


accelerating sclution. 


Perhaps the most significant ground of comparison of 
the two methods other than the shooting problem lies in the 
manner in which each takes advantage of the Hartree- 
Womersley technique. CS merely uses information from _ the 
immediately preceeding x value to obtain initia! 


approximations to the fluid properties: The EJT can use as 
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Many previous x values as necessary to ensure convergence of 
the scheme under certain conditions on the Hessian matrix of 
the system. The EJT makes use of all available information 
to give a certain credibility to the numerical results 
obtained, and a priori information as to the degree of 
truncation error involved. For CS the validity of solutions 
obtained has had to be established by experimental testing. 


2) Programming Considerations. 


The CS is very hard to program and debug due to _ the 
number of complex formulae and the subroutine structures. 
The EJT is simple to program, the only complex formulae 
being involved in the routine setting up the linear 
equations. All other program segments are standard ones, 
for example, one to solve Ryerear of equations. The complex 
routine is made simpler by the use of supporting subroutines 
setting up product and quotient formulae and evaluating 
members cf the Jacobian matrix. The penalty of the easier 
programming is prior manipulation of the problem to put it 
in quasilinearized form but this process is standard and a 
small price to pay. 


3) Computing Time Involved. 


An accurate comparison of the computirg time involved 
3) either method cannot be made without actual 
Weel enereaeton of the EJT . A prediction that the EJT would 
be considerably more efficient can be made, however, based 


on the fact that it is a singly iterative method. 
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The CS results are presented in tabular form at 
iteraticn arguments. Large storage arrays are required for 
this, and the shooting technique involves storage of more 
than one trajectory at each x value until the final version 
is obtained. Supplementary arrays are also needed to test 
for convergence in the iteration loops. The EJT stores a 
set of Chebyshev coefficients at each x value and must 
Maintain a Jacobian matrix at all times. The set of 
coefficients is smaller than the corresponding table of 
values. Interpolation between values in a table for 
arbitrary values of n is avoided. Furthermore, integrals 


or derivatives of stored quantities are readily calculable. 


— 


: 54 a ae. t J [ 
hace’ aa 3 ae 
“id Riandd ale a 


7 
ar) ie . 


ein 3a rte 
soteney (oat? aa paren x one, + See ofonk eat 
fesy oF MefSza os le Sas \sysuze ~ 88 
5B geaots Th2 sd? <agool getyetert ad? ie ak: 
foeu fas onlay Kk doen PS eynelsbitens .¥ 
to.to@ off seewit Usa te “*Erte9 anidorst 
% effet pathneg@etvon oir wit). tofdent =f efm 
a0? fdat 4 ei apiler. aveweid apaastorzeret 
elorgssed | setoutad it — .SsBnees) At.) 30, apy 

sidslppios yitheet 228 aeiedsnenp hervds2 Ic Sovt: 


CHAPTER ELEVEN 


The primary features of the Clutter-Smith (CS) and 
extended Jaffe-Thomas techniques (EJT) techniques will be 


summarized and conclusions drawn. 


The CS method uses the Hartree-Womersley technique to 
transform the equations into ordinary diiferential equations 
in the independent variable n « In that form the equations 
have nenlinear coefficients (the fluid properties) and an 
incomplete initial value set. If the momentum equation is 
to be solved, shooting on the parameter Db" (x, 0) must be 
performed; if this process is successful, the energy 
equation can then be solved as a linear problem in a The 
solution at each value of x is refined by re-calculating the 
coefficients and solving the equations again until a 
convergence criterion is satisfied. A predictor-corrector 
integration technique using four points is used once 
suitable starting formulae have ensured accuracy at the wall 
ie O . The integration formulae can be chosen so that 
error control depends on Lipschitz constants of the solution 
variables, but a priori error control is difficult to ensure 
as these constants are hard to calculate given the 
nonlinearity of the coefficients and the iterative way they 
are treated. Practical values of the fluid parameters in 
the equations make shooting a less than satisfactory 


technique as large variations in trajectory result from 
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} 
Small'@changes “in finitiabiovaluiés lcof D e Various 
improvements on the shooting technique have been proposed 


but the problem cannot be eliminated using the CS approach. 


The EJT method transforms the nonlinear equations into 


@ quasilinearized form which can be integrated using 


techniques for linear equations. The Hartree-Womersley 
technique fits well with the method, as a first 
approximation which ensures convergence of the 


quasilinearization algorithm is provided, as long as 
measures on the soluticn tendencies, which can be estimated 
using the quasilinearization parameters, are within a 
reasonakle range. Thus points of irregularity in the flow 
can be approached more nearly using the EJT than the CS 
method. Advantage can be taken of the linearized form of 
the equation to eliminate the necessity for shooting by 
transforming the differential equations into a set of linear 
algebraic equations by means of a series form of solution. 
Chebyshev polynomials are suitable, and the solution 
representation is attractive in terms of easy maniptulation 
and compact storage. Difficulties encountered an 
representing the dependence of the fluid properties on the 
solution can be overcome uSing interpolation. The method 
solves both the momentum and the energy equations in a 
Single iterative loop for each value of x and the refinement 
of the fluid properties is an automatic consequence of the 


quasilinearization process. The method is easy to implement 
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at the cost of pricribmangpulation to put the problem in 


quasilinear form. 


It must be concluded that the EJT is potentially a 
much mcre powerful technique for solving this kind of 
complex problem than the CS, which is representative of a 
large class of empirically developed methods. Ls 
recommended therefore that the EJT be implemented fully 
using a ‘problem of practical Significance. The a priori 
error ccntrol, avoidance of shooting, and relative ease of 
implementation are incentive enough, while the other 
advantages of the method which have been outlined during the 


thesis are added benefits of the method. 
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The derivation of the equations can be followed in 
detaql ran (105 17,721,720 |. a Prist discussion of the crigan 
of these equations is given below. 

Consider a volume element ina fluid. The mass of fluid 


enclosed in that fluid at time t is 


le spay. 
V 


The total mass change outward due to fluid flow 
operates through the surface of the volume element and is 
; hy Lak 
SE 2) { pU.nds 
where n is the outward normal. The law of conservation of 
mass ina fluid states that the outflow of mass must balance 


the net change of mass inside the element, or 


p ga U.nd 
at hg pay Ib. ih aes 


Using Green's theorem this becomes 
90 4 5 (op i 
i (et i (DU pidVire 0 


and since the volume element is arbitrary the integrand must 


be 0 or 


If the Cartesian coordinates in two-dimensional flow 
Aree teandevpotranstotmatacnhhto thercoordinate esystem ginepr 


described in Chapter One yields 
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Vv. (pu) = 2 (= (rpu) + 2 (rpv)) 


dy 
and making the assumpticn that the flow is steady 
OD 
yoauary 
the continuity equation becomes 
ve 3 3 “ey 
Se (mers + = (rov = 0 
= (er (rpu) ay (rpv) ) 
Zhe momentum equation. 


The Navier-Stokes equation of motion for fluids is derived 


from the expression 


J 
where 
DU _ 3U , = Ux 
Hoe on + U.VU 
Here the second law of Newton is used as follows. The 
integral oi 
DU 
f ae odt 
T 


is the sum of the products of mass and acceleration over the 
volume element T. (The material derivative 2, is 
descriptive of the history of a particular collection of 
particles in the fluid and not of a point in the stationary 
Cartesian space.) The ‘body force per unit mass! or 


"external' force is called F and the integral 


ieee Fdt 
T 
sums these forces over the volume element. It contains al 


external influences acting on the element. Added to the 
effect of these forces is the force of surface tension. A 
stress tensor G.: may be introduced and the sum of the 


t 


effects of these forces may be written 
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Op .n. ass fe: 14 
where the double subscript summation convention is used. 
Newton's law may then be written 


Du; 
{ 5c part = i E,pdt + fs me dt 
T 


Ory, Sintec on abi dar ays 


Substitutions for the material derivative and the 
desired form of the stress tensor yield the momentum 
equation as found in Chapter One. 

The energy equation. 
The energy equation is an expression of the balance of 
thermodynamic energy change in a compressible fluid. It may 
be develcped from the relation 

0 at = go + ae + 3%; (k 3x; 
where ph expresses the energy stored in a volume element 
‘Te Now h 2s enthalpy, defined by 

es le 

where E is the total internal energy of the particles in the 
element. And pt is a term expressing possible variation in 
energy through relative changes in pressure or volume. D 
is the dissipation function or rate of work done by viscous 
forces. ve is the material derivative of pressure with 


respect to time. Also k is the thermal conductivity and T 


the absclute temperature. 


The above equation may be regarded as a statement of 
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the seccnd law of thermodynamics 
des SvdaS y= spay 
in terms of converiently measureable thermodynamic 


quantities. 


The derivation of the specific expression found in 
Chapter Cne from the above form may be followed in Reference 


[28]. 
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TRANSFORMATION OF THE EQUATIONS 


A standard transformation used in boundary layer 

theory 2s that of Dorodnitsvn [19 } 

of 

Tey era oye 

fe) 
which has the effect of stretching the boundary layer 
coordinate normal to the body coordinate proportionally to 
the density, thus making the layer thickness approximately 
uniform along the body. Glutver ~ and "Smith-[-2875. 1] "use »a 
transformation closely related to this one: 
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This transformation iS reminiscent also of the Falkner Skan 
transformaticn f[ 61: 
* Uy el 72 
n (se y 
and allows convenient parallels to be drawn with the 


incompressible equation of momentun. 


Further treatment of the equations consists of 
representing the dependent variables as potential functions 
and thus increasing the order of the equations by one. Such 
treatment is standard in flow studies as it allows automatic 
satisfaction of the continuity eguation by the stream 
function : 
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Clutter and Smith use a dimensionless version of the stream 
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Final transformation of the momentum equation is to 
the quantity 
p= £9 
which helps to reduce the accumulation of round-off error in 
the numerical calculations to be performed. Details of the 
conversion of the momentum equation to its final form may be 


found -une ols 


The energy equation undergoes transformation first to 
the | °{ 1 xX) plane as described above and then to the unknown 
(enthalpy) H. H is transformed by 

g = H/H 
to a dimensionless quantity and then by 
) = g-l 


to cut dcwn round-off error. It is then in its final forn. 
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